import arcpy
from arcpy.sa import *
def InsertAndClip(filepath,insertItem):
    inputfile="C:/Users/xulzu/Desktop/paper/origin_data/csv/"+str(filepath)+"data.csv"
    # 由xy导入数据
    shp=arcpy.MakeXYEventLayer_management(inputfile, "log", "lat", str(filepath)+'data',"", "")
    # 进行插值
    outIDW = Idw(shp, insertItem, 'C:/Users/xulzu/Desktop/paper/origin_data/iamge/imagescale.tif', 2, '')
    # 裁剪后的保存路径
    resultfile="C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\middata\\"+str(filepath)+"data.tif"
    # 按研究区裁剪并保存
    arcpy.Clip_management(outIDW,'#',resultfile,"lanzhoushiwgs84",'',"ClippingGeometry")
for i in range(6):
    inputfile=str(2014+i)
    InsertAndClip(inputfile,'tem')



# ////相关分析
import numpy as np
import arcpy
# 降水栅格数组
rainList=[]
# ndvi栅格数组
ndviList=[]
for x in range(15):
    rainrasterpath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\rainRaster\\'+str(2005+0)+'data.tif'
    rainimage=arcpy.RasterToNumPyArray(rainrasterpath)
    rainList.append(rainimage)
    ndviimagepath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearndvi\\'+str(2005+x)+'ndvi.tif'
    ndviimage=arcpy.RasterToNumPyArray(ndviimagepath)
    ndviList.append(ndviimage)
    
# ////相关分析
# 求影像1和影像2间的相关系数 first 为降水栅格或者温度栅格，second为ndvi栅格
def corrToraster(imageList1,imageList2):
    first=imageList1[0]
    row=first.shape[0]
    column=first.shape[1]
    # 对每一个像素进行操作
    # 初始化一个栅格图
    initRaster=np.zeros( (row,column) )
    for rowindex in range(row):
        for columnindex in range(column):
            # 如果检测到nodata区域则不计算这个位置 

            if(first[rowindex][columnindex]<-200):
                    continue
            fisrtarr=[]
            for image in imageList1:
                fisrtarr.append(image[rowindex,columnindex])
            secondarr=[]
            for image in imageList2:
                secondarr.append(image[rowindex][columnindex])
            x=np.vstack((fisrtarr,secondarr))
            # //得到同一个像元上温度或者降水与ndvi的相关系数阵
            re=np.corrcoef(x)
            # 将计算得到相关系数赋给新栅格
            initRaster[rowindex][columnindex]=re[1][0]
    myRaster = arcpy.NumPyArrayToRaster(initRaster)
    myRaster.save("C:/Users/xulzu/Desktop/paper/origin_data/iamge/rainNdviscale.tif")
corrToraster(rainList,ndviList)

    
# 栅格分割
# C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\temRaster
# C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\temRaster
# C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\rainRaster
# C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\rainRaster
# C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearndvi
# C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\ndviRaster
# 0,1,2,3,4,5,9,10,11,12,13,14,20,21,22,23,30,31,40,68,69,77,78,79,87,88,89,95,96,97,98,99


import arcpy

for i in range(15):
    path='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearsResidual\\'+str(2005+i)+'data.tif'
    basename=str(2005+i)+'data'
    arcpy.SplitRaster_management(path, "C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\15yearclip",
                              basename, "NUMBER_OF_TILES", "TIFF", "NEAREST",
                             "10 10", "#", "#", "PIXELS", "#", "#")


            



